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ABSTRACT 

We present skeleton studies of non-Gaussianity in the Cosmic Microwave Background 
temperature anisotropy observed in the five-year Wilkinson Microwave Anisotropy 
Probe ( WMAP) data. The local skeleton is traced on the 2D sphere by cubic spline 
interpolation which leads to more accurate estimation of the intersection positions 
between the skeleton and the secondary pixels than conventional linear interpolation. 
We demonstrate that the skeleton-based estimator of non-Gaussianity of the local type 
(/nl''') " the departure of the length distribution from the corresponding Gaussian 



yields an unbiased and sufficiently converged likelihood function for 



/NL 

expectation 

flocsd 

/nl ■ 

We analyse the skeleton statistics in the WMAP 5-year combined V- and W-band 
data outside the Galactic base-mask determined from the KQ75 sky-coverage. The 
results are consistent with Gaussian simulations of the the best-fitting cosmological 
model, but deviate from the previous results determined using the WMAP 1-year 
data. We show that it is unlikely that the improved skeleton tracing method, the 
omission of Q-band data, the modification of the foreground-template fitting method 
or the absence of 6 extended regions in the new mask contribute to such a deviation. 
However, the application of the KpO base-mask in data processing does improve the 
consistency with the WMAPl results. 

The /^L^'^'likelihood functions of the data are estimated at 9 different smoothing 
levels. It is unexpected that the best-fit values show positive correlation with the 
smoothing scales. Further investigation argues against a point-source or goodness-of- 
fit explanation but finds that about 30% of either Gaussian or /nl samples having 
better goodness-of-fit than the WMAP 5-year data show a similar correlation. We 
present the estimate /nl**' — 47.3 ± 34.9 (Icr error) determined from the first four 
smoothing angles and /^l^' = 76.8 ± 43.1 for the combination of all nine. The former 
result may be overestimated at the 0.21cr-level because of point sources. 
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1 INTRODUCTION 

Generic inflationary models predict that the initial con- 
ditions of the post-inflation universe can be described 
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by a Gaussian random-phase field witli nearly scale- 
invariant fluctuations. These subsequently seed the per- 
turbations that generate Cosmic Microwave Background 
(CMB) anisotropics and structure formation thereafter. The 
Gaussianity of the statistics determined from measures of 
the CMB anisotropy and large scale structure distribution 
can provide evidence that validates the inflationary scenario 
of the extremely-early Universe. Besides the simplest single- 
scalar fiel d model that predicts a truly Gaussian initial 
condition (lGuthlll98ll: iBardeen. Steinhardt fc TurneijflQS^ : 



iMukhanov. Feldman fc Brandenbergeii Il992l ). there are a 
number of inflationary models predicting non-Gaussianity 
in two broad classifications, the equilateral type and the lo- 
cal type. The detection of a specific type of non-Gaussianity 
can shed light on the fundamental physical properties of in- 
fiation. 

In this paper, we are concerned with a local type non- 
Gaussianity of the "simple st weak nonlinear coupling" case 
JKomatsu fc SpergeJlJOoH ') 
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where $(a;) denotes the primordial curvature perturba- 
tions and $L is its linear Gaussian part. The ampli- 
tude of the non-Gaussianity is parameterised by the di- 
mensionless coupling constant /n^ ' (/nl hereafter). The 
first observational constraint on /nl - —35 00 < /nl < 
2000 at 95 % C.L. - was discussed by iKomatsu et al.l 
l|2002h using the angular bisp ectrum compu ted from the 
four-year COBE DMR data (JBennett et al.| [T996). A re- 
duce d bispectrum technique, hereafter the KSW estima- 
tor, IjKomatsu. Spe rgcl fc Wandclt 2005), was applied to the 
first-year an d three-year WMAP data, leading to —58 < 
/nl < 134 llKomatsu et al.1 12003;1 and -54 < /nl < 114 



l|Spergel et al.ll2007l '). respectively. lYadav fc Wandeltl (120081) 
employed an apparently improved estimator ( Yadav et al.l 
[2008 ) to obtain 27 < /nl < 147 for the V-|-W-band data 
outside the KpO mask with i!max ~ 750 with the three- 
year WMAP data. The WMAP team used the same esti- 
mator to measure /nl from the five-year WMAP V-l-W- 
band outside the KQ75 mask with ^max = 700 and obtained 
-9 < /nl < 111. 

The possibility of detecting CMB non-Gaussianity using 
a group o f morph ological statis tics, - Minkowski function- 
als (M Fs) (JMatsub ara 2003; Hik age. Komatsu fc Matsubaral 
I2OO61 I - has also been studied. The departure of MFs from 
their Gaussian expectations has been tested to be an unbi- 
ased estimator for /nl and then applied to the WMAP3- 
year Q-l-V-l-W combined map yielding —70 < /nl < 91 at 
the 95% C.L. (Hikageet al. 200i). The WMAP team re- 
investigated the MFs estimator with the 5-year template- 
cleaned V-l-W map outside the KQ75 mask, yielding /nl ~ 
-57 ± 60 (68% C.L.) at r esolution N^e = 128 and /nl = 
-68 ± 69 at A^side = 64 IJKomatsu et al.ll2009l ). It is stiU 
unclear why the MFs favour a negative best-fit ampli- 
tude for /nl while the bispectrum estimator prefers a pos- 
itive one, even though the MFs can be formed by the 
weighted sum of the bispectrum. Thus it is of great im- 
portance to use different estimators to identify and inves- 
tigate the weak non-Gaussian signal in WMAP observa- 
tions. In fact, the one-point probability density function 
(1-pdf) of t he smoothed temperatu re field can also be im- 
plemented IJBernardeau et al.l 120021 ) as an alternative non- 



Gaussian ity estimator |jeong fc Smooti 2007^ . Indeed, as 
noted bv lNovikov. Colombi fc Dorg 1 20061 ), the normalised 
differential length of the skeleton is closely linked to this 
quantity, but the skeleton remains of interest due to its dif- 
ferent sensitivity to specific aspects of the data, eg. the noise 
distribution. It is likely that a complete understanding of 
the data can only be realised after the application of a wide 
range of statistical tests. 

The skeleton has been considered as a probe of the fil- 
amentary structures of a 2D or 3D smooth random field. 
The original definition of the skeleton is non-local, making 
the analytical discussion diffi cult and the numerical evalua- 
tion costly. [Novikov, Colom bi fc Dorg (2006) first proposed 
a local approximation that "the local skeleton is given by 
the set of points where the gradient is aligned with the lo- 
cal curvature major axis and where the second component 
of the local curvature is negative". They also presented a 
numerical approach to trace the local skeleton and found an 
approximate expression for the differential length distribu- 
tion of a Gaussian field. As another morphological statisti- 
cal test, the method has been applied to both large-scale 
structure measures ("Sousbi e et"an |2006| . I2OO8I ) and CMB 
anisotropics (Eriksen et al. 20041 ). The latter was performed 
on the Q-l-V-l-W map of the first-year WMAP data out- 
side a base-mask that is defined on the KpO sky-coverage. 
Comparing with Gaussian simulations, the length distribu- 
tion of the skeleton did not show significant deviation from 
the Gaussian predictions. The impact of non-excluded point 
sources was found to be small for the statistics concerned. 

In parallel to studies of non-Gaussian signal estima- 
tors, several algorithms of simulating non-Gaussi an realisa- 
tions have been developed. iKomatsu et al.l ([2003) first sim- 
ulated the local-type non-Gaussian component by integrat- 
ing the spherical harmonics of ^'[(x) — V^^ J d'^x$'l{x) 
in spherical harmonic space. Another strategy has been 
developed in which a pre-computed 'filter' encoding the 
correlation properties of Gaussian curvature perturbation 
multipoles boosts the computation of high-resolution tem- 
perature and pola risation Gaussian and corresponding 
non-Gaussian maps IjLiguori. M atarres e fc Moscardinill2003l : 
iLiguori et al. 20071). This m ethod was recently improved 
bv lElsner fc Wandeltl (|2009l ). Such /nl simulation methods 
provide the community with powerful tools to investigate 
the primordial non-Gaussianity and the impact of other as- 
trophysical and systematic effects on it. 

In this paper, the skeleton length distribution is adopted 
as an estimator of the local-type non-Gaussianity. We adopt 
the cubic spline interpolation to trace the underlying local 
skeleton rather than the conventional linear one to make 
a more accurate estimation of the intersection position be- 
tween the skeleton and pixel edge. Motivated by MFs studies 
on /nl, the statistical properties of the skeleton length dis- 
tribution and the convergence of an /nl estimation method- 
ology are investigated from the /nl simulations. We then 
analyse the skeleton statistics in the five-year release of the 
WMAP data and compare with both Gaussian and non-zero 
/nl samples. The re sults of the nu ll Gaus sian test are com- 
pared with those of lEriksen et al.l l|2004h for the first-year 
WMAP data, and then we use the skeleton estimator to 
compute a likelihood estimate for /nl ■ 

This paper is organised as follows. In Section[2l we carry 
out numerical studies on the CMB local skeleton, includ- 
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ing the skeleton statistics utilised in our analysis (Section 
12. 1[) and the test of unbiasedness and convergency of /nl- 
Ukehhood led by skeleton estimator from noise-free /nl- 
simulations (Section 12. 2|) . [STTI presents an overview of the 
WMAP data and the instrumental properties that should 
be encoded into our simulations to make an unbiased com- 
parison and parameter estimation. Section [3. 21 describes the 
process of computing the estimator and further analysis 
from both the observed data and simulations having con- 
sistent instrumental properties and sky-coverage. Results 
are reported in Section |4l including the analysis and dis- 
cussion of a Gaussian frequentist test (Section l4.ip and /nl- 
estimations (|4.2p . Finally, we present our conclusions in Sec- 
tion [S] 



NUMERICAL STUDIES ON CMB LOCAL 
SKELETON 



According 



to 



the appro ximation made by 



iNovikov. Colombi &: Dorg (|200q ). the local skeleton on 
a smooth 2D sphere p(r), traces those points where the 
gradient of p is the eigenvector of the corresponding Hessian 
matrix. That is, it satisfies the characteristic equation 



nVp = AVp 



(2) 



with A (Ai > A2; A2 < 0) the eigenvalues, where H = 
d'^ p/ drjdrj is the Hessian matrix at position r. Identically 
with lEriksen et al.l (|2004l ). we do not specify the condition 
of eigenvalues of the local linear system. In other words, the 
skeleton in our analysis is considered as the set of underlying 
zero-contour lines of the realisation 



S = PxPyipxx - pyy) + Pxy'ypX 



- Px)-, 



(3) 



where pi and pij denote the first and second derivatives of 
p(r) in two orthogonal directions, x and y. As for the CMB 
temperature field T(n), the 'skeleton map' 5 is re-expressed 
as 



5 = T-eT-(f){T-ee ~ T-, 



4,4,) 



" J-\B4>\J-:<j, — J--,e, 



(4) 



where the semicolons denote the covariant derivatives 
and the definite expressio n of them can be found in 
ISchmalzing fc Gorskil (|2002l ). 

The method for tracing the local skeleton in 
the HEALPix sche me has been reviewed in detail by 
lEriksen et al.l (|2004l ). In Appendix VK[ we seek to optimise 
the method by applying the cubic spline interpolation for 
estimating the underlying positions of skeleton 'knots' on 
the pixelised sphere. The resulting skeleton statistics are in- 
troduced and tested for their applicability to non-Gaussian 
signal detection and /nl estimation. 



2.1 The statistics 

In this work, the CMB temperature realisation intended for 
skeleton analysis, T{n), is first normalised as. 



u{n) = 



Tin) 



(5) 



The standard deviation a is computed over the valid region 
of each realisation after application of an adequate smooth- 
ing process (Section 13. 2|) . 



We utilise the skeleton length distribution function of 
the normalised temperature thresholds v, as a probe of non- 
Gaussianity and to construct an estimator of /nl. As with 
any probability density function, there are two types of dis- 
tributions quantifying the skeleton length, the differential 
pdf 



Cd{u) = 



1 dLjiy) 
I/tot dv 



and the cumulative one 



r + oc 

jC-a{l^) = / Cd{l'')du' , 

J 1/ 



(6) 



(7) 



where the normalisation factor Ltot = f °° dLiv) is the 
total length. 

These two functions are equivalent and should lead to 
consistent results. In the first investigation of the statisti- 
cal properties of th e skeleton length in the WMAP data 
(Erik sen et al.l 120041 ). the cumulative form was utilised and 
compared with the predictions of a Gaussian model. In our 
analysis, both the differential and cumulative functions are 
computed. 

2.2 The idealised skeleton- /wL-test 

We study the signature of the local-type non-Gaussianity as 
a function of /nl on the skeleton length distributions, Cd{i^) 
and Caiy)- As a necessary precursor to /NL-estimation, we 
establish that our estimators lead to an unbiased and suf- 
ficiently converged /NL-likelihood by analysing noise-free 
full-sky realisations with a non-Gaussian signal component. 
The test is based on simulations of the CMB anisotropy 
as a function of /nl . We adopt th e algo rit hm proposed 



by Liguori. Matarrese fc Moscardinil (120031'): Liguori et al.l 



(!2007^ and recentlv improved bv lElsner fc Wandeltl (| 2009| ) 
to simulate a set of Gaussian realisations (a^) with cor- 
responding non-Gaussian components (o^^). The cosmo- 
logical parameters adopted for the /nl simulations are 
those determined for the WMAP5 best-fit ACDM model 
fKomats u et al.ll2009l ) . Specifically, the following parameters 
are adopted: Oa = 0.742, D.ch'^ = 0.1099, Qbh'^ = 0.02273, 
A'^iko = 0.002Mpc-^) = 2.41 x 10"^ h = 0.719, n^ = 
0.963, and r = 0.087. There are a total of 2500 simulated 
{af^, a^} pairs in this test that include power up to a 
maximum multipole fmax ~ 1024. 

Pixelised skymaps with different /nl values are there- 
fore obtained following the relation 

r(p, /nl)=^ ^ (af„ -|-/NLa?m)fe^y«m(p), (8) 

where be is a Gaussian beam transfer function with 
FWHM = 30' and 60' in this test. The first and second 
derivatives of the map can be computed by the HEALPix 
routine almSmap.der. Using the method discussed in Ap- 
pendix [X] the skeleton length distribution /!(:/, /nl) can 
then be estimated from the skeleton map. In this process, the 
normalised temperature threshold is set to ;/ £ [—4.0, 4.0] 
with 25 uniform bins. 

Given the additive nature of the non-Gaussian compo- 
nent, it is reasonable to express Civ) as 



r(z., /nl) = r°M + r^'^C/^, /nl). 



(9) 
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-4.0 -2.0 0.0 2.0 

Threshold level v (ct) 



4.0 



Figure 1. The skeleton estimators obtained from 2500 /nl- 
simulations smoothed by a Gaussian beam of FWHM = 30'. 
Upper: The differential estimator C^ ■ The solid curve depicts 
the averaged distribution, £^'-'(;/, /nl = 150), of 200 bins, with 
the error bars marking the corresponding la errors of the 25 re- 
binned values (black filled circles), while the dashed curve depicts 
the estimation of (£^*^(i^, /nl = —150)). The grey bands corre- 
spond to the 1(t(68.26%) and 2(t(95.44%) confidence regions of 
C^^{u, /nl = 0), i.e., the Gaussian condition. Lower: The cumu- 
lative estimator C^ ■ The nomenclature of elements follows the 
same style as the upper panel. 



For each C{v, /nl) sample, the non-Gaussian component can 
be estimated as 



-C (i^, /nl) 



z:(i.,/nl)-(£^M>, (10) 

where {C {v)) gives the Gaussian expectation of the skele- 
ton length. We depict the samples of C}^ {y, /nl = 0, ±150) 
in Figure [T] The grey bands indicate the la and 2a confi- 
dence regions of a purely Gaussian ensemble, /nl = 0. It 
is noteworthy that the behaviour of the non-Gaussian ex- 
pectation values (£'^ (i/, /nl)) for both the differential and 
cumulative distributions have a characteristic variation with 
threshold. It is similar to MFs in that the peak-trough or- 
der and the amplitude of such features indicate the sign 
and the magnitude of /nl, respectively. This suggests that 
the skeleton can be considered as another morphological 
/NL-estimator, which may lead to deeper understanding of 
the underlying non-Gaussian properties of the observations. 
However, with respect to the la error of C^ , the fluctuation 
is roughly within the 2a range of Gaussian predictions, even 
with /nl = 150 which is larger than the 95% confidence 



The 



best-fit /nl and la error from the likelihood 



2]vEf=iX'(/iLl/NL°)] (Figure EJ computed by differ- 



ential and cumulative estimators from A'^ 



500 noise-free simu- 
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level upper limit for recent /NL-estimations using WMAP 
data. It would still be challenging for a skeleton estimator 
to provide a firm Gaussian/non-Gaussian assessment using 
the observed data. 

Considering only the diagonal elements of the covari- 
ance matrix, we use 2000 simulations to estimate {C {v)), 
the mean and the standard deviation of C^'^{v, fm^). The 
500 remaining simulations are used to compute the x^ func- 
tions. Given a hypothetical value of /nl°) the X^(./nl|/nl'') 
of each /NL-skeleton sample with index i (i = 1, 2, ...500) is 
computed as 



i ( ri I rtruc 

UnlI/nl 



)-E 



cn^^h^) -{£-'''' {v,r^r)) 



<7(£NG(^J^.uc)) 



(11) 

where the correlations between bins have not been taken 
into account because the full covariance matrix is not suf- 
ficiently converged for the available sample volume in our 
analysis. Further tests indicate that the corresponding like- 
lihood from each sample is of bimodal or even multi-modal 
shape if the full covariance matrix is adopted, which causes 
the estimation to be unrevealing. 

The parameter /nl is uniformly sampled from —300 to 
300 with a step-length A/nl = 5. We estimate the likeli- 
hoods for three specific /nl'^ values, and ±150. The pos- 
terior PDF for /nl° can be obtained by Bayes' theorem 

p(/Nri{,fNL}) oc p({/^L}i./sr) X p{f^r) 



Wnn 



I /'truc\ 

nlI/nl ) 



(12) 



oc exp 



'o^Z^^^-^nlL/nD 



where we have conservatively set the prior P{flii^°) to be 
uniform and A'' equals to 500. In fact, we have found that 
roughly 20 samples with FWHM= 30' smoothing are ad- 
equate for the posterior distribution to converge sharply 
around /nl"^ with a la error ~ A/nl = 5. However, we have 
only one observed CMB sample so that the convergence of 
the consequent posterior distribution is limited by the data 
resolution and the noise level. The effective likelihood func- 
tions of each sample, i.e., exp [^-^^ J2'iLi X^(/nl1/nl'')J , are 
illustrated in Figure [21 using different normalisation factors 
for visual convenience. The histograms depict the computed 
likelihoods which are perfectly fitted by Gaussian functions. 
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Figure 2. cxp [-^ E£i X^(/nlI/nl")] . the effective likeli- 
hood functions computed by the differential estimator C^ with 
input parameter f^^" = 0, ±150. The actual functions are renor- 
malised by different factors for visual convenience and are shown 
by histograms. They are extremely well fitted by the Gaussian 
functions depicted by solid curves. The higher, narrower (lower, 
wider) histograms and curves correspond to likelihoods from sim- 
ulations with FWHM = 30' (60'). 



Accordingly, the mean and the la width of each likehhood 
are estimated as presented in Table [l] 

The results demonstrate a good recovery of the input 
/nl values given the interval A/nl = 5 of our sampling. The 
X (/nlI/nl'^) in Eg. 1111 therefore constitutes an unbiased 
maximum likelihood position in /NL-space and the corre- 
sponding la error is determined by the likelihood function. 
It is noteworthy that the cumulative estimator behaves a 
little bit better than the differential one and therefore the 
former is selected for /nl estimation as applied to real data. 



3 METHOD 

Even though the literature contains theoretical predictions 
for the length distributions of the local skeleton on a 2D 
Gaussian random field, our analysis compares measures de- 
rived from simulated observations of the sky with the corre- 
sponding values for the WMAP data, since the inhomoge- 
neous noise contribution and the complicated sky-coverage 
render analytical investigation difficult. Furthermore, it is 
also difficult to interpret the non-Gaussian component of 
these skeleton measures analytically. In what follows, we 
introduce both the instrumental properties impacting the 
observed data and the essential numerical processing steps 
required for further analysis. 



3.1 The WMAP data and the simulations 

The WMAP instrument measures the CMB temperature 
anisotropy in five fre quency bands from 23 to 94 GHz 
IjBennett et al.l l2003al ) . The foreground-reduced sky maps 
in V and W-band are used in our analysis, identical to the 
data sel ection for the WM AP five-year power spectrum esti- 
mation (jNolta et al.ll2009l ). These maps are available in the 
HEALPix pixelisation scheme with A'aido ~ 512 from the 




t.^_^i.- 






■^r\_^^ 




Figure 3. Comparison of the KQ75 base-mask (KQ75B) used 
in WMAP5 data processing with the KpO base-mask (KQ75B) 
used for WMAPl and WMAP3 analysis. The white (light-grey) 
regions are excluded (included) by both KQ75B and KpOB. The 
dark-grey (black) parts are excluded by the former (latter) but 
not excluded by the latter (former). There are 6 extended regions 
(labeled from '1' to '6') eliminated in the KQ75B mask. 

LAMBDA websitqj. The maps from two (four) differencing 
assemblies (DAs) at V- (W-) band are combined using uni- 
form weights over the sky and equal weights for each DA. 
The resulting maps in V- and W-band are then combined 
to obtain the VW-band map using the same method. The 
effective beam transfer function of the VW-band can then 
be easily computed from the beam functions of those DA^ 
constituting the VW-band map. The observational data are 
inevitably affected by the instrumental noise, dominated by 
an uncorrelated component with a variance per pixel de- 
pending on the noise amplitud e ap and the pixel sc anning 
strategies of each DA, Naha{p) (jBennett et al.ll2003al ). 

The extended temperature analysis mask (KQ75) is 
adopted to minimise the contamination from the diffuse 
Galactic foreground and point source emission. For further 
investigations, the part related to the Galactic emission is 
separated out to form a base-mask called 'KQ75B' in our 
analysis. As a comparison of the base-masks used in the 
1-, 3- and 5-year WMAP data analyses, we illustrate the 
KpOB (adopted by tEriksen et al. (J004)) and KQ75B mask 
in Figure O Besides the extended Galactic profile, there are 
six extended regions (labeled from '1' to '6') eliminated by 
the new base-mask. The impact of these regions on skele- 
ton statistics will be considered w hen comparing our results 
with those of lEriksen et al.l l|2004l ) for WMAPl. 

The 2500 pairs of Gaussian and non-Gaussian realisa- 
tions {afrnjO-^m} mtroduccd in Section [2.21 are used for our 
/nl studies. For each /nl value, we construct a map with 
resolution parameter A'^sidc ~ 512 and WMAP instrumental 
properties as 



T{p, /nl) 



y^ y^ {afm + fNhafm)bePeYem{p) 

e=2 m=-t 

aoV 



iVcbs(p)' 



(13) 



where bi is the effective beam transfer function of the 

t http://lambda.gsfc.nasa.gov/product/map/dr3/ 

maps_da_forered_r9_iqu_5yr_get . cfm 

^ http://lambda.gsfc.nasa.gov/product/map/dr3/beam_info.cfm 



© 2010 RAS, MNRAS OOO.ITHTtI 



6 Hou et al. 



WMAP VW-band data and pe is the pixelisation window 
function for A'aide = 512. The second term on the rhs sim- 
ulates the noise contribution on each pixel with Gaussian 
random number rj ~ A'^(0, 1). 

In this work, we perform both a Gaussian frequentist 
test and /NL-estimations. In the former, the Gaussian sim- 
ulations are processed in the same way as Eq.[T3]but free of 
the /nl term. 

3.2 Data processing and the analysis 

In this section, we introduce the data processing methods 
applied to both the observed and the simulated realisations 
for studies of the skeleton length distribution. The process- 
ing steps presented here follo w the strategy detailed in Sec- 
tion 4 of lEriksen et al.1 (|2004 ). 



3.2.1 Map processing 

The base-mask is applied to the map to avoid Galactic 
foreground contami nation. Following the methodology of 
lEriksen et al.l (12004! ) ■ we do not exclude point sources, in 
particular because any additional smoothing applied to the 
mask reduces the sky coverage available for analysis dramat- 
ically. This approach is supported by studies of the spec- 
tral parameter, 7, bv iEriksen et al.l ([2004), which indicates 
that smoothing of the data renders the skeleton less sensi- 
tive to point source signal contributions for larger FWHMs. 
Moreover, a median-filter technique is applied to the point 
sources to investigate their impact on the skeleton statis- 
tics for smaller smoothing FWHMs. Specifically, for a given 
pixel i that would be eliminated by the point-source mask, 
we consider all other unmasked pixels within a 1° radius and 
determine the median temperature for this set of pixels. The 
temperature at pixel i is then replaced by this median value, 
and the process repeated for all pixels specified in the point- 
source mask. The median-filtered map is then analysed in 
the same manner as the unfiltered data set. 

Following standard procedure in CMB data analysis, 
the monopole and dipole components are fitted and re- 
moved from each map outside the masked region. This 
step is achieved using the HEALPix F90 subroutine, re- 
move.dipole. The resulting map is then smoothed with a 
Gaussian beam in harmonic space, again using HEALPix 
tools. For easy com parison with the previous analysis of 
lEriksen et al] (|2004h . the FWHM widths, 6Ifwhm, selected 
for these Gaussian smoothing beams are taken to be 0?53, 
0?64, 0?85, 1?28, 1?70, 2?13, 2?55, 2?98 and 3? 40 (we aban- 
don the larger angular scales used in the former analysis 
to ensure good convergence of the /NL-likelihood). Since a 
higher resolution map is necessary for more accurate esti- 
mation of the skeleton statistics (see Appendix IX|) . the res- 
olution parameter of the resulting smoothed map is set to 
Af.ido = 1024. 

After the smoothing of the data, the base-mask must 
also be expanded and the same processing method is fol- 
lowed. For each smoothing scale, only those base-mask pix- 
els with values larger than 0.99 are defined to be valid pixels 
on the smoothed mask. 

Finally, on the valid region defined by the smoothed 
base-mask, the map T{p), either from observation or simu- 
lation, is renormalised to temperature thresholds, v (Eq[5]), 



while the invalid pixels are abandoned for computing the 
standard deviation. Using the method discussed in Ap- 
pendix [X] the skeleton length distributions, Cd{v) and 
Caiy), can be estimated for each set of smoothed samples. 
The original distribution Liu) in Eq. |6] is divided into 200 
bins with v G [—4.0, 4.0] during skeleton tracing. 

3.2.2 Non-Gausstan detector and estimator 

From the processed Gaussian simulations, we compute the 
Gaussian expectation of the skeleton statistics for each 
smoothing scale, {£,'^{v, 6'pwhm)). The departure from these 
expectation values is then obtained for both the observed 
data and each Gaussian sample as 



A£,{u,9fwiim) ~ C{iy, Ofwhm) 



{C°{u, 



0), (14) 



and the corresponding x^ value is then computed 



2 
X = 



A£,{u,6fwhm) 



_a{A£.{u,eFWHM))_ ^^^^ 

where we omit the (AC) term since it is definitely 
zero. In the /nl analysis, the non-Gaussian departure, 
£ (:/, /nl, fewHM), and the X statistics are estimated by 
Ea llOl and llll The best-fit value and error of /nl can then be 
obtained by analysing the likelihood function as discussed 
in Section [2^21 

Before we provide final estimates of /nl from the differ- 
ent smoothing scales, we combined the estimators, AjC{v) of 
the data and C {v, /nl) of each set of /nl sample, to 



ACciiyjNL) 



JVtwhm 

y~^ Wi(i^, /nl)A£'(i/) 



and 



Affwh,, 



-Co (J'Jnl) 



^ Wi{u,ftii,)£^'~^'\v,fN 



respectively with the inverse-variance weighting 



Wi{v,fNL) 



Erir"i/'^^(^^'^''(^,/NL)) 



(16) 



(17) 



(18) 



where i corresponds to one smoothing scale and Aifwhm rep- 
resents the number of scales used in the combination. The 
combined x is then computed 



Xc(/n 



^f A£c(/^,/NL)-{£g'^(^,/NL)) , .^„. 

4^1 a[£NG(^^^^^)] I ■ i^y) 



This combination makes an integrated estimation of /nl 
which includes the non-Gaussian signal at several different 
scales with a mild weighting. 



4 RESULTS AND DISCUSSIONS 

4.1 Gaussian frequentist results 

We first compare the observed results with our Gaussian 
model predictions. In this case, we perform 10240 Gaussian 
simulations of the WMAP VW-band properties. Different 
base-masks, as well as the median-filter, are applied indepen- 
dently to both the real and the simulated skies to study the 
foreground effect on the skeleton results. The corresponding 
X^ values are then computed to enable the frequentist test. 
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4.1.1 Results of KQ75B processing 

For each smoothing scale, the skeleton length departure from 
the Gaussian expectation, AjC{v, Sfwhm) = jC{v, Sfwhm) — 
{j-PiVjO-pwiiM)), is computed from samples obtained with 
the KQ75B masked maps. The results are shown in the left 
two columns (for both the differential and cumulative dis- 
tributions) of Figure E] for Sfwhm = 0?64, 0?85, 1?28, 1?70, 
2?98 and 3? 40. The grey bands demonstrate the la and 2a 
confidence regions of the Gaussian prediction. The observed 
ones are rebinned to 25 bins and depicted by filled circles 
with the Icr-error bar of each bin. The rebinning is necessary 
since the differential skeleton distribution is relatively noisy. 

In the case of the cumulative distributions, A£,a{j^) for 
WMAP5, some features consistent with a positive /nl value 
are observed, albeit within the la Gaussian confidence level. 
The behaviour of the differential distribution, A£tj(!/), sup- 
ports this inference despite the existence of a higher level of 
fluctuations. However, there are differences be tween the new 
results and the corresponding WMAPl ones (lEriksen et al.l 
[2OOJ). For each smoothing scale, the latter show a Icr- level 
peak around u = while the neighbouring troughs show less 
fluctuations especially in the 1/ > 1 region. In contrast, as 
shown in Figure [5] (the left two columns), the former's peak 
is less apparent but the troughs are much more distinct par- 
ticularly for 6'fwhm = 1?28 and 1?70. The comparison be- 
tween WMAPl and our new results is shown in Figure [7] for 
6'fwhm = 0?64, 0?85 and 1?28. 

There are several possibilities associated with such a 
discrepancy. 

(1) Change of the skeleton-tracing method. Utilising cu- 
bic spline interpolation in the skeleton tracing algorithm 
yields a more accurate estimation of the quantities than the 
previously adopted linear algorithm (see Appendix |B]| . We 
computed /S.Ca{v) for the template-cleaned WMAPl data 
us ing the same band-se lection, mask and processing steps as 
in lEriksen et al.l l|2004l ). and tracing the underlying skeleton 
by both linear and cubic spline interpolation strategies. The 
Gaussian expectation is also estimated in both cases using 
simulations. The results are shown in the left column of Fig- 
ure [7| where it can be seen that the dashed black line (cubic 
spline) and the solid grey line (linear) essentially overlap. 
It is therefore clear that changing the interpolation scheme 
contributes little to the discrepancy found. This issue is also 
discussed in Appendix [Bl for the WMAP5 data; 

(2) Band-selection. In the analysis of lEriksen et al.l 
(1200J), the Q-, V- and W-band maps are combined with a 
spatially-invariant inverse-noise-variance weigthing. The re- 
sulting map is dominated by the Q-band since it has the 
lowest noise of the three. However, since it is the band 
for which Galactic foreground residuals remain significant, 
it is plausible that these have an impact on the skeleton 
results. We repeated our analysis using the appropriately 
weighted WMAP5 Q-, V- and W-band data, but retain- 
ing the KQ75B base-mask. Corresponding Gaussian simula- 
tions are also performed. The results are shown as the black 
connected-filled-circles in the right column of Figure [7] The 
profile shows modest deviation from our VW-results (black 
fiUed-squares), however, it does not result in the discrepancy 
level required. On the contrary, the difference becomes less 
significant for large Sfwhm. 

(3) Difference of the foreground subtraction method be- 




t. -u^x^^fiWSSi- 



0.030 [mK] 



Figure 4. The difference map between the WMAPl and 
WMAP5 combined foreground models at V-band. between 
WMAPl and WMAP5. The KpOB and KQ75B masks are also 
denoted. The ringing effect around some bright sources, especially 
the L MC, comes from the o ne-year processing of the three tem- 
plates JBennett et al.ll2003bl 'l. 



tween WMAPl and WMAP5. The foreg round templates 
used for the former (iBennett et al.ll2003bl 'i are the FDS 94 
GHz dust prediction, the Ha map for free-free emission 
and the 408 MHz Haslam map for synchr otron emission. 
The three-year WMAP foreground analysis IjHinshaw et al.l 
120071 ) and beyond replace the 408 MHz data with a tem- 
plate based on the the K-Ka difference map. The differ- 
ence between the two foreground models at V-band utilis- 
ing the coefficients for the first-year fits of Bennett et al.l 
(|2003bt ) and the five-year analysis of ICold et all (|2009l ') is 
shown in Figure |4] The profile demonstrates a dipole-like 
structure in the large-scale temperature distribution out- 
side both the KpOB or KQ75B masks, which may affect the 
skeleton statistics and the corresponding inferences of /nl. 
We subtract the five-year foreground model from the one- 
year raw maps at Q-, V- and W-bands , which are then com - 
bined and processed identically with lEriksen et al.l (|2004r ) 
using the KpOB mask. The corresponding skeleton statis- 
tic, A£,a{i^), is depicted by the connected open-circles in 
the left column of Figure [T] They demonstrate consistency 
with the original WMAPl results. Similarly, another inde- 
pendent test has been carried out on the five-year raw maps 
from which the one-year foreground model is subtracted be- 
fore the data are combined and processed using the KQ75B 
mask. The results are depicted as the dashed grey line in the 
nght column of Figure [T] and demonstrate consistency with 
our five-year templated-cleaned VW-KQ75B results (black 
fiUed-squares). We conclude that it is difficult to attribute 
the observed discrepancy to the change of foreground sub- 
traction method. 

(4) Change of the base-mask in processing. It is very sus- 
picious that the residual foreground components around the 
dark-grey regions in Figure [3] bias the skeleton results of 
WMAPl, although mild smoothing and mask thresholding 
are applied before skeleton tracing. We discuss this issue 
in Section 14.1.21 by investigating the Galactic plane region 
and the extragalactic sources (labeled from 1-6 in Figure [Sjl 
separately. 



As shown by the left two columns in Figure [5] the pro- 
file of the WMAP5 AjC-a{v) function is consistent with that 
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Table 2. The processing elements for the styles of lines and sym- 
bols in Figure [71 



Band 


Mask" 


Fore- red'' 


Intcrp'^ 


Style 


QVW 


KpOB 


lyr 


Cubic 


Dashed blacl< line 


QVW 


KpOB 


Syr 


Cubic 


Connected open-eircles 


QVW 


KQ75B 


lyr 


Cubic 


Dashed grey line 


QVW 


KQ75B 


Syr 


Cubic 


Connected filled-cireles 


QVW 


KpOB 


lyr 


Linear 


Solid grey line 


VW* 


KQ75B 


Syr 


Cubic 


Solid black line 



^ The base-mask applied in map-processing and analysis. 

^ The templates and the corresponding coefficients applied for 
foreground-reducing before our map-processing. 

■^ The interpolation method used for tracing the underlying 
local skeleton. 

'^ The V-l-W combined data with uniform weighting, while spa- 
tial invariant inverse-noise-variance weighting for QVW. 



Table 3. The x^-based frequentist results for the WMAP5 skele- 
ton analysis derived using different processing masks and meth- 
ods on 10 smoothing scales. We list the fraction of the simulations 
with a x^ values less extreme than the observed one. The letters 
'M' correspond to 'median-filter'. The values are determined from 
the estimator ACa{i^) computed for 200 bins from the data and 
10240 Gaussian sam ples. The corresponding results for WMAPl 
JEriksen et al.ll2004l ) are also listed for easy comparison. 



FWHM WMAPl KQ75B KQ75M KpOB KQhybrid 



0?53 


0.234 


0.1220 


0.1115 


0.0812 


0.1310 


0?64 


0.286 


0.1503 


0.1345 


0.1539 


0.1604 


0?85 


0.354 


0.2608 


0.2148 


0.2147 


0.2720 


1?28 


0.293 


0.3490 


0.3167 


0.2481 


0.3590 


1?70 


0.284 


0.4258 


0.3761 


0.1360 


0.4504 


2? 13 


0.248 


0.3691 


0.3745 


0.1669 


0.3822 


2? 55 


0.208 


0.3205 


0.3352 


0.1379 


0.3361 


2? 98 


0.166 


0.2728 


0.2684 


0.1343 


0.2892 


3? 40 


0.113 


0.2119 


0.2389 


0.1023 


0.2237 


3? 83 


0.081 


0.1866 


0.2410 


0.0923 


0.1963 



expected for a positive- /nl- In particular, both ACd and 
A£a, rebinned for FWHM = 1?28 and 1?70, demonstrate 
consistent features with the solid lines shown in Figure [l] 
for /nl = +150. However, the troughs in the v > 1 region 
(hot region) seem relatively less depressed. It is likely that 
the point sources and foreground components contribute to 
this asymmetry between the two troughs. The results from 
the median-filtered map yield insights implications into this 
issue. 

We computed the x^ values of A£a for both the ob- 
served and the simulated samples. We list the fraction of 
the simulations with a x values less extreme than the ob- 
served one in Table [S] Th e corresponding WM APl results 
are also listed (Table 3 in lEriksen et al.l (|2004l ')'). Generally 
speaking, there is no qualitative difference between the five- 
year and one-year results. But our results show a unimodal 
dependence on the smoothing scales. The /NL-signal seems 
more significant around the angular scales FWHM = 1?28, 
1?70 and 2? 13. 



unchanged. We also create a new base-mask called 'KQhy- 
brid' which excludes the same Galactic plane with KQ75B 
but handles the six extended sources (Figure (H]) identically 
to KpOB. The KQhybrid mask is then included in the data 
processing too as an independent test. Some of the results 
are shown in Figure [5] and the right column of Figure [T] 

In general, the ACa profiles of the KpOB processing are 
generally consistent with the previous WMAPl results, al- 
though the peak-trough structure is not identical in detail. 
The KQhybrid mask yields a consistent set of results with 
those of KQ75B as shown in Figure [B] Moreover, we have 
applied the KQ75B mask to the WMAPl data and found 
that the results (dashed grey line in the left column of Fig- 
ure [T} show a similar discrepancy from the KpOB processed 
ones and consistency with results from our 5-year VW data 
processing. This indicates that modifications of the mask do 
significantly affect the skeleton estimation in the WMAPl 
analysis. Although the reason can be easily found by exam- 
ining the area ratio of the dark-grey regions in Figure [3] it is 
important to make a separate investigation on the impact of 
residual Galactic foreground and extragalactic sources since 
/nl analysis exhibits different res ponses to diffe rent types of 
foreground contamination t Cabella et al.ll2010r ). This sepa- 
rate analysis motivates future skeleton studies on the effects 
of different Galactic foreground templates. 

It is noteworthy that the skeleton discrepancies caused 
by base-mask selection indicate that residual Galactic fore- 
grounds bias the non-Gaussian analyses for WMAPl and 
even WMAP3 since the KpO mask was the standard temper- 
ature analysis window then and the Kp2 mask excluded even 
less area around Galactic plane. This issue may have impli- 
cations on the bispectrum analysis because the additional 
smoothing operation, which smears the local structures of 
foreground templates, is not necessary for bispectrum esti- 
mation. 

The foreground issue is also assessed as a comple- 
ment to the mask-changing analysis. We subtracted the five- 
year (one-year) foreground templates from the raw maps of 
WMAPl {WMAP5) data. The subtracted maps are then 
combined and processed using the KQ75B (KpOB) mask 
and the skeleton results are depicted as the connected fiUed- 
circles (dashed black line) in the left (right) column of Fig- 
ure [7] They are consistent with the results from the stan- 
dard foreground subtraction processing with the same corre- 
sponding base-mask. It is therefore confirmed that the fore- 
ground model is not responsible for the discrepancy of the 
skeleton statistics as seen. 

The corresponding results are listed in Table |3] It is 
straightforward to infer that the KQhybrid processing re- 
sults are more consistent with the corresponding KQ75B 
ones. The differences of a few percent come from the 6 ex- 
tended regions. The x^ results from the KpOB-processing are 
somewhat different to the WMAPl inference although the 
profiles of Ajda are quite similar. Besides the band selection, 
it is most probably due to the modified template-fitting of 
the Galactic foreground in five-year data processing, as well 
as the better S/N level in 5-year data. 



4-1.2 Results of KpOB and KQhybrid processing 

We ap plied the one- year KpOB mask used in lEriksen et al.l 
(I2OOJ) in our analysis with all other operations remaining 



4.1.3 Results of median-filter processing 

In order to assess the validity of the results from the KQ75B 
mask analysis in which point sources are not excluded, we 
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KQ75base 



KQ75 median-filter 




KpObase 



KQtiybrid 



-2.0 0.0 2.0 4.0 -2.0 0,0 2.0 4.0 -2,0 0.0 

Threshold level (v) in a 

Figure 5. The skeleton statistics, A.C{u) = C{u) — {C{v)), com- 
puted from KQ75B (left two columns) and KQ75 median-filter 
(right two columns) processing on smoothing scales FWHM = 
0?64, 0?85, 1?28, 1?70, 2?98 and 3? 40. The first and third column 
(the second and forth column) correspond to results of the differ- 
ential (cumulative) estimator. The grey bands show the la and 2a 
confidence regions defined by 10240 Gaussian samples. The black 
fiUed-circles connected by solid lines show the observed sample 
which is rebinned to 25 bins. The error bar mark the la error in 
each bin according to such rebinning. 



have applied a median filter to those pixels located at po- 
sitions in the point source mask before smoothing, then 
processed the filtered map to obtain the skeleton statistics. 
Some results are plotted in the right two columns of Figure 
O and listed in Table [3] In general, the median filtered re- 
sults show good consistency with the KQ75B results even 
for the first few smoothing scales, implying that the base- 
mask processing is safe for skeleton analysis on the scales 
considered in this work. 

Nevertheless, small visual differences suggest further 
investigation into how point sources modify the skeleton 
statistics and the /nl estimations. We make a comparison 
of the WMAP?) ACa between the KQ75B and median-fiher 
processing. The differences between them are plotted in Fig- 
ure [8] for FWHM = 0?53, 0?64, 0?85 and 1?28 by solid, dot- 
dashed, dashed and dotted-lines, respectively. It is suggested 
that the point sources do have asymmetric impacts on the 
skeleton for positive and negative temperature thresholds - 
negative biasing is seen for the range —2.0 < v < 0.0 and 
positive biasing is apparent for 0.0 < i^ < 2.0. In particular 




-4.0 -2.0 0.0 2.0 4.0 -2.0 0.0 2.0 4.0 -2,0 0,0 2,0 4,0 -2.0 0.0 2.0 4.0 

Threshold level {v) in a 

Figure 6. The skeleton statistics, A.C{u) = C{u) — {C{u)), com- 
puted from KpOB (left two columns) and KQhybrid (right two 
columns) processing on smoothing scales FWHM = 0?64, 0?85, 
1?28, 1?70, 2?98 and 3? 40. The nomenclature of the elements fol- 
lows the same style as Figure [5] 



for the dotted- line, a 30% lower depression is observed over 
0.0 < I' < 0.5. This could bias the best-fit /nl value though 
the bins around this range are assigned lower weights ac- 
cording in the combined y^ computation. Although the plot 
suggests that the magnitude of potential biasing seems to 
increase with smoothing scale, the larger smoothing still re- 
duces sensitivity to point sources. Moreover, the profiles seen 
in Figure |8] become increasingly noise-like within the range 
I' G [—2.5, 2.5] at larger smoothing scales. 



4.2 /n 



estimation 



4-2.1 General results 

Using the method introduced in Section [2.21 the likelihood 
function for /nl is estimated for each smoothing scale based 
on the 2500 sets of /nl samples, £'^'~'(i/, {/nl}). We sam- 
ple the parameter within the range /nl € [—200,400] with 
step-length A/nl = 2.5. The KQ75B-processed data are 
utilised from FWHM = 0?53 to 3? 40 with the median-fiher- 
processed data from FWHM = 0?53 to 1?28 compared for 
reference. We use the cumulative estimator A£a because it 
leads to 10% more converged estimations than the differen- 
tial one according to a mock test (Section 12. 2p . Before x^ 
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WMAPl 



WMAP5 




2.0 4.0 -2.0 

Threshold level (i/) in a 



Figure 7. The skeleton statistics, ACa{i^), computed from diflFerent processing methods. Left Column: the results obtained from the 
one-year WMAP data. Right Column: the results obtained from the five-year WMAP data. The lines and symbols denoting different 
processing elements are noted in Table [2] 



computation, the estimator resulting from both the observed 
data and simulations are rebinned to 25 bino 

The results are shown in the top panel of Figure [9] with 
each curve depicting the likelihood (without normalisation) 
for each smoothing scale. The likelihood functions are fitted 
by Gaussian functions so that the best-fitting /nl and the 
corresponding la error are obtained and then marked in the 
same plot. The likelihood at the highest resolution indicates 
that the Gaussian hypothesis (/nl = 0) is rejected only at 
0.8cr-level, while it increases to 2.7cr for FWHM = 2? 13. It 
is apparent that the best-fitting /nl values show a positive 
correlation with the smoothing scale, which is unexpected 
since /nl is scale-independent according to the local-type 
non-Gaussian model and our simulations. 

As discussed in Section [4.1.31 although the estimation 
is inevitably biased by the point sources or other types of 
foreground, large angle smoothing renders the estimation 
insensitive to those effects. We repeat the estimation us- 



■^ It has been tested that 25 is the best number for rebinning 
in our analysis. More bins will make the estimator more noisy 
so that the resulting likelihood is bimodal or even multimodal, 
whereas less bins will make the likelihood less converged. 



ing median-filtered samples from the first four smoothing 
scales. As shown in the middle panel of Figure [H] the re- 
sults are consistent in general, and the positive correlation 
between /Nif* and the smoothing scales is identical to the 
unfiltered analysis. It is therefore suggested that the point 
sources contribute little to such correlation. The la errors 
are robust according to the median-filter reference but the 
best-fit values of /nl from the KQ75B processing seem to 
be over-estimated by levels of 0.04(t, 0.26(t, 0.39(t and 0.22(t 
for FWHM = 0?53, 0?64, 0?85 and 1?28, respectively 

In principle, different heights of the the likelihoods rep- 
resent variations in the goodness-of-fit if the correspond- 
ing x^ values have the same number of degrees-of-freedom. 
A higher likelihood implies the /nl expectation fits the 
data better and it does appear that the likelihoods from 
larger-angle smoothing (FWHM = 2?55, 2?98 and 3? 40) 
show better results than for smaller FWHMs. However, 
in our analysis, we pick up only the diagonal elements of 
the covariance matrix to compute the x^- It is inappropri- 
ate to make a theoretical interpretation of the goodness- 
of-fit. Consequently, the correlation found above would be 
a false appearance because there might be some bad fit- 
tings. For each FWHM, the x^ value at the maximum likeli- 
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-2024 
Threshold level {v) in a 

Figure 8. The distribution difference between KQ75B and KQ75 
median-filter processed estimator, A£a of WMAP5 data. The 
cases of different smoothing scales are distinguished by different 
line-styles. 



Table 4. The goodness of fit, i.e., the probabilities that the sim- 
ulated samples with {xL„(/NLl/Si^)} < X^in- f^l^ and xLn 
are the maximum likelihood /nl of the data and its correspond- 
ing x^ value of each case, respectively. The results in the case of 
combined data (KQ75B Comb.) are also listed. 



FWHM 



KQ75B 



P(%) 



KQ75M 

din Pi%) 



KQ75B Comb. 



0?53 


6.50 


27.1 


6.28 


20.9 


0?64 


6.34 


25.8 


6.57 


22.3 


0?85 


7.89 


35.6 


8.18 


30.5 


1?28 


7.85 


33.5 


7.88 


27.3 


1?70 


9.42 


40.5 


N/A 


N/A 


2? 13 


6.15 


19.5 


N/A 


N/A 


2° 55 


5.20 


12. 1'' 


N/A 


N/A 


2° 98 


5.26 


11.7 


N/A 


N/A 


3? 40 


4.37 


6.0 


N/A 


N/A 



;.80 



;.24 



30.4 



28.0 



' This number may be under-estimated because the underlying 
X^ minima of some samples lay outside our /nl sampling range, 
i.e., their corresponding /j^j^ > 400. Similar cases are also found 



for FWHM = 2? 98 and 3? 40. 



hood (ML) of the data is represented as Xmin- Accordingly, 
there are 2500 sampled X^({/nl}|/nl') and each has a min- 
imum within our sampling range. We count the probabil- 
ity of xiin({/NL}|/NL') < X^in to quantify the goodness of 
fit for results from both the KQ75B and median-filter pro- 
cessing, with a lower probability corresponding to a better 
fit. The Xmin values and the probabilities are listed in Ta- 
ble m The moderate probabilities are consistent with each 
other though they may be under-estimated for the last three 
FWHMs. On one hand, it is demonstrated that our skeleton 
statistic fits the possible /nl feature in the WMAP5 data 
and our estimations are therefore validated. On the other 
hand, it remains unconfirmed what the source of the posi- 
tive correlation between /niT's and smoothing scales is, and 
which we will return to in Section 14.2.31 



i. 



i. 



u.ia 




KWHM (deg) 


AT' I" 


. . 1 , 




' 








0.53 


25.7 30.9 






/ 


\ " 


0.10 




0.64 

0.85 

1-28 


37-5 33-1 
65-3 37-8 
112 49-9 








0,08 


- 


1-70 

2-13 


156 65-4 
230 85-4 








....^ " 


0.06 


- 


2-55 

2-98 

3-40 


274 103 
302 121 
304 138 


/ 


"^ \~ 


0.04 






/' \^ 




A 




■■, 


0.02 






A /A 


xy-^ 


\ 




4-^ 


0.00 




^; 



-200 



-100 



100 200 

/nl 



300 



400 



0-OB 














FWHM (deg) 
0-53 


i5r' 'i" ' ' ' '\ 

24.4 31-0 : 


0-05 








/- 


■> 




0-64 

0-85 

1-28 


28-7 33-0 - 
50-5 37-8 : 
101 49-8 ; 


0,0-1 


- 






/ 




\ 




- 








1 


/ 




\ 




: 


0-03 


; 




1 1 
1 1 








\ 


1 


0-02 

0-01 
00 


- 


..^ 


1 1 

1 1 
t 
li 
// 




y 


y 




\ : 
\ - 



-100 



is 



0-02 r 



0-01 




200 



-200 -100 



300 



Figure 9. The likelihood functions of /nl from the skeleton 
statistic for WMAP5 data. Top: The /nl likelihood functions 
computed by KQ75B processed ACai^) and C^ {i^, /nl) oil 9 
different smoothing scales. The estimator is rebinned to 25 bins 
before analysis. The best-fittings and la errors are obtained by 
fitting the likelihoods using Gaussian functions. Middle: Sim- 
ilar cases for statistic derived by KQ75 median-filter process- 
ing on 4 smoothing scales. Bottom: The /nl likelihood func- 
tions estimated from the combined estimator A£a,c(!^i /ng) and 
'^o' (^i/ng)- The solid (dot-dashed) curve shows the likelihood 
computed from KQ75B processed estimator with 4 (9) FWHMs 
combined. The dashed curve corresponds to the KQ75 median- 
filter processed likelihood. 
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4-2.2 Estimation from the com,hined /S.Ca 



As presented in Section 13.2.21 the combinations on differ- 
ent smoothing scales are apphed separately to the rebinned 
/S.Ca{i') of the data and ^^'^(i^, /nl) of the /nl samples. It 
is verified that such a combination still leads to an unbiased 
estimation of /nl (Appendix [C|. 

In our analysis, the first 4 and all 9 scales are combined, 
yielding estimates of /nl,c = 47.3 ±34.9 and /nl,c ~ 76.8 ± 
43.1 respectively, by fitting the likelihood using a Gaussian 
function. The likelihoods are shown in the bottom panel 
of Figure |9] and the goodness-of-fit is also listed in Table 
13] The estimates are consistent with the results discussed 
in Section [4. 2. II and the moderate probabilities (30.4% and 
28.0%) validate the best-fit results. 

The median-filtered results are also combined over the 
first 4 FWHMs and the corresponding likelihood is depicted 
by the dashed curve, resulting in the estimate /nl.c = 
39.8 ± 34.9. The point sources lead to an over-estimate of 
/niT.c ^t the 0.21(T-level according to this comparison. The 
combined estimators, A£a,c(«^) for the KQ75B processed 
data and C^%(i', /nl ~ 0, 47.5, 77.5Jj for the corresponding 
/nl simulations, are illustrated in Figure [Tn]for comparison. 



4-. 2. 3 Cosmic variance and /niT* 

It is interesting that /Nif*s shows a monotonic correlation 
with smoothing scale. The discussions above argue against 
the explanation based on point sources or goodness-of-fit. 
We search for this kind of correlation in our mock samples 
to investigate whether cosmic variance is a possible source 
of such a correlation. In order to make a comprehensive and 
reliable interpretation, we pick up those Gaussian and /nl 
samples which show /nl features at least to the same extent 
as the WMAP5 data. The selection method is introduced 
below. 

(1) Gaussian samples. Similar to the WMAP5 data, each 
of the 10240 Gaussian samples of ACa{i^), is input into /nl- 
estimations on all 9 FWHMs as introduced in Section [32121 
The chi-square for each FWHM, XGauss(/NL,^FWHM), is ob- 
tained as a function of smoothing scale and /nl before we 
combine the estimators of all 9 FWHMs to AjCcii^, Inl)- 
The combined chi-square, Xc,Gauss(/NL), and likelihood are 
then computed by the combined estimator. We find 3111 
samples whose min[xc, Gauss (/nl)] are less than Xc.min from 
the WMAP5 data. It is believed that these samples demon- 
strate better /NL-like features (of /JJl') than the WMAP5 
data for all 9 smoothing scales even though there is no non- 
Gaussian component encoded in the simulations. 

(2) /nl samples. For the 9- FWHM combination discussed 
in Section|3221 samples with Xc,mi„(/NL|/^^L"° = /Sl.c) less 
than the WMAP5 Xc min 8.re selected from 2500 groups of 
/nl samples. The 701 selected samples form the /nl ref- 
erence for investigating the correlation between /Sif^s and 
smoothing scale. 

In 3111 selected Gaussian samples, we find 844 that 
feature a monotonic correlation with smoothing scale (~ 



Note that the step-length for /nl sampling is 2.5 in our anal- 
ysis, /nl = 47.5 and 77.5 are the maximum likelihood values. 
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Figure 10. The combined skeleton estimator of the data, Gaus- 
sian and non-Gaussian predictions with AfpwHM = 4 (upper) 
and 9 (lower). The thin-solid curves show ACa.c('^) of WMAP5 
data with the filled-circles showing the corresponding 25 bins' re- 
binnings. The grey bands represent the la and 2cr confidence re- 
gions of the Gaussian predictions of ACa ci'^)- According to 2500 
/nl simulations, the thick-solid curves depict the expectation of 
£^^(1^, /nl = 47.5,77.5) for Afp-yvHM = 4 and 9, respectively, 
and the dashed and dot-dashed curves depict the corresponding 
la and 2a confidence boundaries. 



27.1%). Similarly, there are 222 /nl samples from 701 show- 
ing the same behaviour ('^ 31.7%). According to our tests, 
they show similar properties to that illustrated in the top 
panel of Figure [9] where the maximum likelihood for large- 
scale smoothing is 'pulled' significantly to the non-Gaussian 
region. There is a considerable probability (around 30%) of 
such a correlation so that cosmic variance is a highly prob- 
able explanation. 



5 CONCLUSIONS 

In this paper, we have studied the local-approximation to the 
skeleton on a 2D sphere pixelised in the HEALPix scheme, 
and refined the method of tracing the quantity. The statis- 
tical properties of the skeleton estimator have subsequently 
been investigated using mock CMB temperature anisotropy 
maps. 

The cubic spline interpolation method locates the skele- 
ton knots more accurately than the simple linear method, 
which makes the local linear system more robust at the 
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knots. This is of great importance for finer analysis of the 
local system. For example, the studies on skeleton classifica- 
tion (|Pogosvan et al.ll2009l ). which is performed by analysing 
the eigenvalues of the linear characteristic equation, request 
highly accurate estimation of such eigenvalues in particular 
around the demarcation point between two types of skeleton. 
Our modification provides a more reliable basis for this kind 
of study. The departure of the skeleton length distribution 
from its Gaussian expectation shows connections with both 
the sign and the magnitude of /nl so that it would yield 
a /NL-likelihood function. Based on simulated sets of CMB 
temperature anisotropy with a local type of non-Gaussian 
component, it has been tested that both the differential and 
cumulative skeleton estimators provide unbiased and suffi- 
ciently converged likelihood function for /nl, but the latter 
yields a likelihood 10% more converged than the former. 

The estimator was applied to the five-year WMAP data 
release and the results compared with not only the Gaussian 
predictions, but al so the results f r om th e first-year WMAP 
data processing bv lEriksen et al.l l|2004 ). An /NL-likelihood 
function has been estimated by computing the x^ on the 
basis of 2500 sets of /nl samples. We have also investigated 
the goodness of fit, the impact of the point sources and the 
comic variance effect on the best-fit amplitudes of /nl- The 
analysis is carried out on the V+W combined map for vari- 
ous sky coverages. 

The process i ng ste ps in our analysis follow closely those 
of lEriksen et al.l (|200J) but utilise the new five-year KQ75 
mask and combined V- and W-band data. The smoothing 
scales adopted i n our data processin g are also identical to 
those selected in lEriksen et al.l (1200J). Our skeleton results 
show an apparent deviation from the first-year ones. Accord- 
ing to an extensive series of tests, it is the difference between 
the two Galactic plane regions defined by the KQ75 and KpO 
masks that contributes mostly to the shifts. Generally, the 
KQ75 mask excludes a more extended region close to the 
Galactic plane than the KpO mask, and this should be more 
conservative for temperature analysis. This kind of devia- 
tion to the skeleton estimates implies a systematic bias in 
/NL-estimation, in other words, previous /nl studies carried 
out on KpO sky coverage (or even the less conservative KP2 
mask) may be biased by the residual Galactic foreground 
within the dark-grey regions as shown in Figure [S] We do 
not exclude the pixels located in point sources to allow suffi- 
cient convergence of the likelihood. The impact of the point 
sources on the estimator is investigated by analysing the dif- 
ference to samples using median-filtered maps. The results 
show that the point sources do have an asymmetric impact 
on the estimator between the positive and negative tempera- 
ture thresholds on the four smallest smoothing scales. How- 
ever, the effect is less significant for larger FWHMs. The 
results of our frequentist analysis show that the WMAP5 
data are consistent with Gaussian predictions. 

We have estimated the /nl likelihoods on 9 smoothing 
levels. The results show an unexpected positive-correlation 
between the best-fit amplitudes, /nl*) ^'id FWHM smooth- 
ing scales. The peak of the likelihood function seems to be 
'pulled' to a highly non-Gaussian region with the Gaussian 
case, /nl ~ 0, being 'expelled' to the very tail of likelihood 
for some large smoothing scales. Further investigations argue 
against a point source explanation since the median-filtered 



data still exhibit such a correlation. However, the presence 
of point sources may yield an over-estimation of /n?*- 

The combination of samples for the first 4 and all 9 
smoothing scales lead to the best-fit amplitudes with la 
errors, /nl = 47.3±34.9 and /nl ~ 76.8±43.1, respectively. 
The median-filter studies suggest that the best-fit over 4 
scales may be over-estimated at the 0.21(T-level because of 
point sources. An investigation has been carried out on the 
unexpected correlation between /Sif^s and smoothing scales 
using both Gaussian and /nl samples with a goodness-of-fit 
better than that for WMAP5. About 30% of them show the 
behaviour seen in our analysis, so that cosmic variance may 
be an appropriate explanation for this issue. 
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APPENDIX A: THE LOCAL SKELETON IN 
HEALPIX FRAME 

We construct a local coordinate system on the 2D HEALPix 
sky map shown in Figure lAll where the direction to 
the Galactic north-pole is depicted as 'N'. Following the 
HEALPix coordinate conventionqj, two orthogonal axes, x 
and y in Eq. [31 are set to be aligned with the polar-angle 6 

and azimuth (j) axes, respectively. 

Our starting point is identical with that of lEriksen et al.l 

I1200J) in that we determine a pair of vertices on the edge 
of the pre-constructed secondary pixels on <S, with one 
vertex value lower but the other higher than zero (can- 
celing vertices, hereafter). It is suggested that the under- 
lying skeleton crosses over those edges connecting pairs 
of canceling vertices. Figure lAll illustrates an exagger- 
ated version of this process. Interpolation is then neces- 
sary to determine the positions of the intersections on edges 
(skeleton knots, her eafter). Line ar interpolati on has been 
adopted previously ( Novikov. Fe ldman &: Sha ndarinI 1 19991 : 
IShandarin et al.ll2002l : lEriksen et al.ll2004i '). since it has been 
widely employed in morphological studies on both the CMB 
and large scale structures, eg. the length and genus quan- 
tities of MFs, which are related to the contour lines of the 
random fluctuation field and its derivatives. However, the 



see 'The HEALPix Primer' in software package, version 2.10 



accuracy of an interpolation method is limited by the topo- 
logical properties of the random field and the pixel size of the 
corresponding realisation. Linear interpolation is accurate 
enough for an analysis of the MFs at current observational 
resolutions (A'^sidc ~ 512, 1024), however, it is inadequate to 
provide precise positions of 'knots' on the skeleton map in 
Eq.|4]as a higher-order (cubic) random field. This may intro- 
duce not only bias in to the statistics of the skeleton length 
for a specific realisation, but also could result in a false de- 
termination of the eigenvalue of the local linear system, in 
particular around the demarcation point between the two 
types of skele tons considered f or studies of skeleton classi- 
fications (Pog osvan et al.ll2009l ). Given the cubic nature of 
the skeleton field, we therefore apply a cubic spline inter- 
polation in this analysis, as introduced in the following text 
in detail. A comparison between the linear and cubic spline 
strategy is presented in Appendix [Bl 

Once a pair of canceling vertices has been found, e.g., 
X2 and X3 in Figure lAll the 6 pixels are then picked up 
with canceling vertices in the middle, as xo, xi, ..., x^. The 
connection lines of the 6-pixel centres must cross over the 
pairs of opposite sides of the quadrangular pixels and be 
parallel with the connection line of canceling vertices (e.g., 
12X3). The values of these 6 pixels (vertices of secondary 
pixels), yi = S{xi){i — 0, 1, 2, 3, 4, 5), are utilised to find the 
spline functions along the connection lines, 



S{x) = { 



' So{x) X e [a::o,a;i] 
Si{x) X £ [xi,a;2] 



(Al) 



,5*4(3;) X £ [x4,a;5] 



where each Si is the piecewise cubic polynomial between the 
pixel-centres 



S,{x) 



Zi+l{x - Xi)^ + Zi{Xi+i - x)^ 

6hi 
+ (1^- yz. I {x,+i-x} 



(A2) 



hi is equal to \xi+i — Xi\ corresponding to the radial dis- 
tance of the two pixel-centres. The coefficients {zi} can be 
obtained by solving the linear system 

hi-iZi-i+2{hi-i + hi)zi + hiZi+i — 

20 = Z5 = 0. 

where yi corresponds to the pixel value at Xi, i.e., the skele- 
ton value, S{xi), in this work. 

Note that this 6-point system on the sphere has been 
approximated by a ID straight line since the pixel-size in 
our analysis is so small (A^'sidc ~ 1024, Spix ~ 3.44'). In fact, 
we only need 5*2 (x) to determine the locations of the knots, 
e.g., p3 in Figure [XT] by solving the cubic equation 



S2ix) = 0. 



(A4) 



There is one and only one real root, Xk, satisfying the con- 
dition X2 < Xk < X3. Then the vector of the underlying knot 
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Figure Al. Tracing the underlying contour (dot-dashed curve) on the grid of the HEALPix pixelisation scheme. The Xi mark the 
positions of the primary pixels utilised in the cubic spline interpolation (Eg. lAl l and lA2t . The thick solid lines connect the estimated 
intersection points (pis) between the underlying contour and the edge of the secondary pixels by means of a cubic spline interpolation. 
The definition of secondary pixels is identical to that of lEriksen et all 1120041) . We select a coordinate system consistent with the HEALPix 
convention, which is specified by a North direction, 'N', and two orthogonal directions (d,cj>) for derivative. 



can be obtained as 



Xk 



X3 - Xk Xk - X2 
X2 H X3, 



X3 - X2 



X3 — X2 



(A5) 



and the corresponding temperature value at x^ is 



n = Ei^liT2 + ^J^^^n. 



X3 — X2 



X3 - X2 



(A6) 



According to Figure lAll after determining the vector of p2 
and P3 (i.e., Xk2 and x^a), the skeleton length within the 
secondary pixel, X2X3x'3x'2, can be estimated by the dot- 
product of these two vectors. 



SL{Ts 



Xk2 

\Xk2\ 



Xk3 
\Xk3\ 



(A7) 



The corresponding temperature value of this piece of skele- 
ton length, Ts, is approximately the simple average of Tk2 
and Tk3- 

It is always the case that the four edges of one secondary 
pixel are connecting canceling vertices. Most of these cases 
indicate a stationary point (maxima or minima or saddle 
point) within this secondary pixel, implying two skeletons 
cross inside. There are still a few exceptions but they will be- 
come very rare due to the small pixel-size and the smoothing 
applied afterwards. We th erefore make the same assumption 
as in Erikscn ct al.| (1200J) that all of the cases indicate a pair 
of skeletons crossing over each other. The possible deviation 
from the length distribution is totally negligible according 
to various tests. 



APPENDIX B: COMPARISON BETWEEN 
LINEAR AND CUBIC SPLINE 
INTERPOLATION FOR SKELETON ANALYSIS 

On a pixelised 2D random field, the key step in tracing the 
local skeleton is to locate the skeleton knot which is always 
within the line connecting the centres of the two canceling 
neighbouring pixels (one edge of the secondary pixel) , and 
whose position is conventionally estimated by linear interpo- 
lation, since the skeleton realisation 5 can be considered as a 
linear function along the line connecting just a few pixels at 
a very high resolution-level. This is an approximation that 
makes things easier to handle, especially for the HEALPix 
pixelization scheme. However, the skeleton is actually a cu- 
bic function, so that it is necessary to test whether linear 
interpolation is sufficient for its computation. In this ap- 
pendix, we investigate the linear properties at the skeleton 
knots derived by linear and cubic spline interpolation meth- 
ods. 

The characteristic equation (Eq. [2]) for the 2D random 
field must be satisfied at the skeleton knots. It can be reex- 
pressed for a CMB temperature field as 






We define 



ri 



We ' 






r-2 



A 



T..t 



(Bl) 



T;j,gT-a -j- T-SthT-a 



and A should satisfy the following 



T-ee — A 



T-Bif, 

T-^cj, — A 







(B2) 



with two real roots Ai and A2 (Ai ^ A2). In principle, ri 
should be equal to r2 and also equal to Ai or A2 along the 
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Figure Bl. The pieces of profile of 5 witliin a six-pixel-array re- 
produced by linear and cubic spline interpolation from one Gaus- 
sian simulation. The centres of the six pixels are identified by xq, 
xi, ..., X5 where the values of <S are marked by filled triangles. 
The positions of the skeleton knots estimated by cubic splines 
(linear lines) are located by the filled (open) circles. The small 
plots inside zoom the curves within X2 and X4. The pixel index 
(position) of Xi is selected to be the same for FWHM = 30' and 
60'. 
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Figure B2. The difference of ACa of WMAPb data between the 
cubic spline (cub) and linear (lin) interpolation processing. The 
cases of different smoothing scales are distinguished by different 
line-styles. 



underlying skeleton. However, in practice, we have to investi- 
gate such a property at the position of the estimated skeleton 
knots on the pixelised sphere where r\ and r2 are not ex- 
actly equal because of estimation errors. The first and second 
derivatives there can be obtained safely by linear interpola- 
tion as in Eg. lAGI We define a new quantity r = (ri +r2)/2. 
The numerical robustness of the equivalence between r and 
A indicates the quality of the estimation method. 

In this test, we pick up two six-pixel-arrays 
(a;o, a;i, ..., 15) from one simulated Gaussian realisation (res- 
olution parameter A'sidc ~ 1024) smoothed by Gaussian 
beams with FWHM = 30' and 60'. The pixel location of the 
two arrays are exactly the same with each other. The corre- 
sponding values of S in Eq. 3] are marked by filled triangles 
in Figure iBl] For the case of FWHM — 30', X2X3 is a pair 
of canceling pixels and Pf point (P point) is the estimated 
skeleton knot determined by a linear (cubic spline) interpo- 
lation method. The linear properties at the two points are 
quantified as 



Pf- n = 


-0.6887, 


^2 = 


-0.7466, r = -0.7174 


Ai = 


-0.6422, 


A2 = 


-0.7221 


P: ri = 


-0.7099, 


r-2 = 


-0.7116, r = -0.7108 


Ai = 


-0.6125, 


A2 = 


— 0.T107 (secondary skeleton) 



It is shown in this example that the cubic spline interpola- 
tion leads to a more accurate location of the skeleton knots, 
and the distribution of skeleton length therein. Note that 
there are two suspicious skeleton knots within X3 and X4 in 
this case but they would not be involved in analysis since X3 
and X4 are not canceling pixels. It is also noteworthy that 
the point P belongs to a piece of the first-type sec ondary 
s kelet on according to the classification in Pogosva n et al.l 
(|2009^ . The robust equivalence between r and the eigenvalue 
indicates accurate and unbiased classification, in particular 
around the underlying demarcation point between two types 
of skeleton where the two eigenvalues are quite close to each 
other. The cases for FWHM = 60' are listed below 

P/i : n = 0.0428, ra = -0.4198, r = -0.2313 

Al = -0.2114, A2 = -0.8062 
Pi : n = -0.1916, r2 = -0.2425, r = -0.2170 

Al = -0.2108, A2 = -0.8060 (primary skeleton) 



/2 



n = 


-0.9264, 


r2 = 


-0.7597, r = -0.8430 


Ai = 


-0.2230, 


A2 = 


-0.8091 


n = 


-0.7981, 


^2 = 


-0.8004, r = -0.7993 


Ai = 


-0.2588, 


A2 = 


—0.7996 (secondary skeleton) 



For the difference between the five-year and one-year 
skeleton processing, we must investigate the impact of 
method selection on the results. Given the KQ75B processed 
data and Gaussian simulations, we carry out the skeleton 
analysis following the steps described in Section[3]2]but util- 
ising linear interpolation to locate the skeleton knots. The 
resulting length departure of the data is then obtained 



'-^'-a,WMAP5 — •'-a,Vl 



(r°'"' 



(B3) 



and the differences between the cubic spline and linear re- 
sults are plotted in Figure |B21 for FWHM = 0?53, 0?64, 0?85 
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Table CI. The Maximum-Likelihood (/Jfj^), best-fitting (/^^^f*) 
values and la error from the likelihood ^cC/nlI/J^l") (Figure 
ICl|l computed from the combined estimator derived from A^ = 
250 KQ75B processed noisy simulations with given parameter 



A^FWHM 
4 


ftruc 

/nl 
0.0 


fMh 
JNh 

-2.5 


fbcst 
JNL 
-2.0 


35.5 


9 


0.0 


-2.5 


-1.7 


42.3 
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200.0 


200.0 


199.4 


36.5 


9 


200.0 


197.5 


197.8 


43.4 



Figure CI. exp | -^ Ef^i X^W^lI/nl" ) 



the effective like- 



lihood functions computed by the combined accumulative esti- 
mator £^Q from KQ75B processed noisy simulations with input 
parameter /J^l"^ = 0, 200. The functions are renormalised by dif- 
ferent factors for visual convenience and are shown by histograms. 
The Gaussian-fitting functions are depicted by solid curves. The 
higher, narrower (lower, wider) histograms and curves correspond 
to the combinations with A''pwHM = 4 (9). 



and 1?28. It is noteworthy that the magnitude of such a dif- 
ference contributes less than 10% to the discrepancy between 
the WMAP5 and WMAPl skeleton length distribution pro- 
file. However, the structure shown in Figure lB2l suggests that 
the linear method would lead to an over-enhanced peak and 
over-depressed trough, which for the positive- /nl structure 
of ACa suggested by the data may bias the best-fitting value 
of /nl. 



APPENDIX C: TEST OF THE LIKELIHOODS 
FROM THE COMBINED ESTIMATOR 

In this section we test for the presence of bias in our com- 
bined estimator. Given simulated noisy realisations from 
the KQ75B processing and the predetermined expectation 
(£c (i/, /nl)), we randomly pick up TV = 250 sets of /nl- 
samples, £g°(zv,/^L) U = 1,2,. ..,250) with A^fwhm = 4 
and 9, to form the conditional x^ functions 



Xc(/nlI/nl°) — / , 



^^""{r^Ji^J 



(^g^'(^,/sr)) 



'^[^g°(^,/i^r)] 



(Cl) 



and the effective likelihood function for each sample. 



■^c(/nl|/nl°) ocexp 



Weplot^c(/NL|/,^r) 



1 2 ■ 

T7 / ,Xc(/nlI/i 



2N 



truc^ 
NL J 



J = l 



(C2) 



values (0 and 200) in Figure lCTl for A^'pwhm = 4 and 9, notic- 
ing that the sampling width A/nl is 2.5. Again, the likeli- 
hoods are perfectly fitted by Gaussian functions with the pa- 
rameters listed in Table [CT] Despite the noise contribution 
and sky-cut, it is demonstrated that the inverse-variance- 
combination still leads to an unbiased skeleton estimator 
for /nl- 



© 2010 RAS, MNRAS OOO.ITHTtI 



